home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / liboctave / dbleGEPBAL.cc < prev    next >
C/C++ Source or Header  |  1996-03-03  |  6KB  |  230 lines
  1. /*
  2.  
  3. Copyright (C) 1996 John W. Eaton
  4.  
  5. This file is part of Octave.
  6.  
  7. Octave is free software; you can redistribute it and/or modify it
  8. under the terms of the GNU General Public License as published by the
  9. Free Software Foundation; either version 2, or (at your option) any
  10. later version.
  11.  
  12. Octave is distributed in the hope that it will be useful, but WITHOUT
  13. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  14. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  15. for more details.
  16.  
  17. You should have received a copy of the GNU General Public License
  18. along with Octave; see the file COPYING.  If not, write to the Free
  19. Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  20.  
  21. */
  22.  
  23. #if defined (__GNUG__)
  24. #pragma implementation
  25. #endif
  26.  
  27. #ifdef HAVE_CONFIG_H
  28. #include <config.h>
  29. #endif
  30.  
  31. #include <cmath>
  32.  
  33. #include <string>
  34.  
  35. #include "dbleGEPBAL.h"
  36. #include "f77-fcn.h"
  37.  
  38. extern "C"
  39. {
  40.   int F77_FCN (dgebak, DGEBAK) (const char*, const char*, const int&,
  41.                 const int&, const int&, double*,
  42.                 const int&, double*, const int&, int&,
  43.                 long, long);
  44.  
  45.   int F77_FCN (reduce, REDUCE) (const int&, const int&, double*,
  46.                 const int&, double*, int&, int&,
  47.                 double*, double*);
  48.  
  49.   int F77_FCN (scaleg, SCALEG) (const int&, const int&, double*,
  50.                 const int&, double*, const int&,
  51.                 const int&, double*, double*,
  52.                 double*);
  53.  
  54.   int F77_FCN (gradeq, GRADEQ) (const int&, const int&, double*,
  55.                 const int&, double*, int&, int&,
  56.                 double*, double*);
  57. }
  58.  
  59. int
  60. GEPBALANCE::init (const Matrix& a, const Matrix& b, const string& balance_job)
  61. {
  62.   int a_nr = a.rows ();
  63.   int a_nc = a.cols ();
  64.   int b_nr = b.rows ();
  65.  
  66.   if (a_nr != a_nc || a_nr != b_nr || b_nr != b.cols ())
  67.     {
  68.       (*current_liboctave_error_handler)
  69.     ("GEPBALANCE requires square matrices of the same size");
  70.       return -1;
  71.     }
  72.  
  73.   int n = a_nc;
  74.  
  75.   int info;
  76.   int ilo;
  77.   int ihi;
  78.  
  79.   Array<double> cscale (n);
  80.   double *pcscale = cscale.fortran_vec ();
  81.  
  82.   Matrix wk (n, 6, 0.0);
  83.   double *pwk = wk.fortran_vec ();
  84.  
  85.   // Back out the permutations:
  86.   //
  87.   // cscale contains the exponents of the column scaling factors in its 
  88.   // ilo through ihi locations and the reducing column permutations in 
  89.   // its first ilo-1 and its ihi+1 through n locations.
  90.   //
  91.   // cperm contains the column permutations applied in grading the a and b 
  92.   // submatrices in its ilo through ihi locations.
  93.   //
  94.   // wk contains the exponents of the row scaling factors in its ilo 
  95.   // through ihi locations, the reducing row permutations in its first 
  96.   // ilo-1 and its ihi+1 through n locations, and the row permutations
  97.   // applied in grading the a and b submatrices in its n+ilo through 
  98.   // n+ihi locations.
  99.   
  100.   // Copy matrices into local structure.
  101.  
  102.   balanced_a_mat = a;
  103.   balanced_b_mat = b;
  104.  
  105.   double *p_balanced_a_mat = balanced_a_mat.fortran_vec ();
  106.   double *p_balanced_b_mat = balanced_b_mat.fortran_vec ();
  107.  
  108.   // Check for permutation option.
  109.  
  110.   char job = balance_job[0];
  111.  
  112.   if (job == 'P' || job == 'B')
  113.     {
  114.       F77_XFCN (reduce, REDUCE, (n, n, p_balanced_a_mat, n,
  115.                  p_balanced_b_mat, ilo, ihi,
  116.                  pcscale, pwk));
  117.     }
  118.   else
  119.     {
  120.       // Set up for scaling later.
  121.  
  122.       ilo = 1;
  123.       ihi = n;
  124.     }
  125.  
  126.   if (f77_exception_encountered)
  127.     (*current_liboctave_error_handler) ("unrecoverable error in reduce");
  128.   else
  129.     {
  130.       Array<double> cperm (n);
  131.       double *pcperm = cperm.fortran_vec ();
  132.  
  133.       // Check for scaling option.
  134.  
  135.       if ((job == 'S' || job == 'B') && ilo != ihi)
  136.     {
  137.       F77_XFCN (scaleg, SCALEG, (n, n, p_balanced_a_mat, n,
  138.                      p_balanced_b_mat, ilo, ihi,
  139.                      pcscale, pcperm, pwk));
  140.     }
  141.       else
  142.     {
  143.       // Set scaling data to 0's.
  144.  
  145.       for (int i = ilo-1; i < ihi; i++)
  146.         {
  147.           cscale.elem (i) = 0.0;
  148.           wk.elem (i, 0) = 0.0;
  149.         }
  150.     }
  151.  
  152.       if (f77_exception_encountered)
  153.     (*current_liboctave_error_handler) ("unrecoverable error in scaleg");
  154.       else
  155.     {
  156.       // Scaleg returns exponents, not values, so...
  157.  
  158.       for (int i = ilo-1; i < ihi; i++)
  159.         {
  160.           cscale.elem (i) = pow (2.0, cscale.elem (i));
  161.           wk.elem (i, 0) = pow (2.0, -wk.elem (i, 0));
  162.         }
  163.  
  164.       // Initialize balancing matrices to identity.
  165.  
  166.       left_balancing_mat = Matrix (n, n, 0.0);
  167.       for (int i = 0; i < n; i++)
  168.         left_balancing_mat (i, i) = 1.0;
  169.  
  170.       right_balancing_mat = left_balancing_mat;
  171.  
  172.       double *p_right_balancing_mat = right_balancing_mat.fortran_vec ();
  173.       double *p_left_balancing_mat = left_balancing_mat.fortran_vec ();
  174.  
  175.       // Column permutations/scaling.
  176.  
  177.       char side = 'R';
  178.  
  179.       F77_XFCN (dgebak, DGEBAK, (&job, &side, n, ilo, ihi, pcscale,
  180.                      n, p_right_balancing_mat, n, info,
  181.                      1L, 1L));
  182.     
  183.       if (f77_exception_encountered)
  184.         (*current_liboctave_error_handler)
  185.           ("unrecoverable error in dgebak");
  186.       else
  187.         {
  188.           // Row permutations/scaling.
  189.  
  190.           side = 'L';
  191.  
  192.           F77_XFCN (dgebak, DGEBAK, (&job, &side, n, ilo, ihi, pwk,
  193.                      n, p_left_balancing_mat, n,
  194.                      info, 1L, 1L));
  195.  
  196. #if 0
  197.           // XXX FIXME XXX --- these four lines need to be added and
  198.           // debugged.  GEPBALANCE::init will work without them, though, so
  199.           // here they are.
  200.  
  201.           if ((job == 'P' || job == 'B') && ilo != ihi)
  202.         {
  203.           F77_XFCN (gradeq, GRADEQ, (n, n, p_balanced_a_mat, n,
  204.                          p_balanced_b_mat, ilo, ihi,
  205.                          pcperm, pwk));
  206.         }
  207. #endif
  208.  
  209.           if (f77_exception_encountered)
  210.         (*current_liboctave_error_handler)
  211.           ("unrecoverable error in dgebak");
  212.           else
  213.         {
  214.           // Transpose for aa = cc*a*dd convention...
  215.  
  216.           left_balancing_mat = left_balancing_mat.transpose ();
  217.         }
  218.         }
  219.     }
  220.     }
  221.  
  222.   return info;
  223. }
  224.  
  225. /*
  226. ;;; Local Variables: ***
  227. ;;; mode: C++ ***
  228. ;;; End: ***
  229. */
  230.